home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DMCHB.FOR < prev    next >
Text File  |  1985-11-29  |  10KB  |  295 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DMCHB
  5. C
  6. C        PURPOSE
  7. C           FOR A GIVEN POSITIVE-DEFINITE M BY M MATRIX A WITH SYMMETRIC
  8. C           BAND STRUCTURE AND - IF NECESSARY - A GIVEN GENERAL M BY N
  9. C           MATRIX R, THE FOLLOWING CALCULATIONS (DEPENDENT ON THE
  10. C           VALUE OF THE DECISION PARAMETER IOP) ARE PERFORMED
  11. C           (1) MATRIX A IS FACTORIZED (IF IOP IS NOT NEGATIVE), THAT
  12. C               MEANS BAND MATRIX TU WITH UPPER CODIAGONALS ONLY IS
  13. C               GENERATED ON THE LOCATIONS OF A SUCH THAT
  14. C               TRANSPOSE(TU)*TU=A.
  15. C           (2) MATRIX R IS MULTIPLIED ON THE LEFT BY INVERSE(TU)
  16. C               AND/OR INVERSE(TRANSPOSE(TU)) AND THE RESULT IS STORED
  17. C               IN THE LOCATIONS OF R.
  18. C           THIS SUBROUTINE ESPECIALLY CAN BE USED TO SOLVE THE SYSTEM
  19. C           OF SIMULTANEOUS LINEAR EQUATIONS A*X=R WITH POSITIVE-
  20. C           DEFINITE COEFFICIENT MATRIX A OF SYMMETRIC BAND STRUCTURE.
  21. C
  22. C        USAGE
  23. C           CALL DMCHB (R,A,M,N,MUD,IOP,EPS,IER)
  24. C
  25. C        DESCRIPTION OF PARAMETERS
  26. C           R      - INPUT IN CASES IOP=-3,-2,-1,1,2,3  DOUBLE PRECISION
  27. C                          M BY N RIGHT HAND SIDE MATRIX,
  28. C                          IN CASE IOP=0  IRRELEVANT.
  29. C                    OUTPUT IN CASES IOP=1,-1  INVERSE(A)*R,
  30. C                           IN CASES IOP=2,-2  INVERSE(TU)*R,
  31. C                           IN CASES IOP=3,-3  INVERSE(TRANSPOSE(TU))*R,
  32. C                           IN CASE  IOP=0     UNCHANGED.
  33. C           A      - INPUT IN CASES IOP=0,1,2,3  DOUBLE PRECISION M BY M
  34. C                          POSITIVE-DEFINITE COEFFICIENT MATRIX OF
  35. C                          SYMMETRIC BAND STRUCTURE STORED IN
  36. C                          COMPRESSED FORM (SEE REMARKS),
  37. C                          IN CASES IOP=-1,-2,-3 DOUBLE PRECISION M BY M
  38. C                          BAND MATRIX TU WITH UPPER CODIAGONALS ONLY,
  39. C                          STORED IN COMPRESSED FORM (SEE REMARKS).
  40. C                    OUTPUT IN ALL CASES  BAND MATRIX TU WITH UPPER
  41. C                           CODIAGONALS ONLY, STORED IN COMPRESSED FORM
  42. C                           (THAT MEANS UNCHANGED IF IOP=-1,-2,-3).
  43. C           M      - INPUT VALUE SPECIFYING THE NUMBER OF ROWS AND
  44. C                    COLUMNS OF A AND THE NUMBER OF ROWS OF R.
  45. C           N      - INPUT VALUE SPECIFYING THE NUMBER OF COLUMNS OF R
  46. C                    (IRRELEVANT IN CASE IOP=0).
  47. C           MUD    - INPUT VALUE SPECIFYING THE NUMBER OF UPPER
  48. C                    CODIAGONALS OF A.
  49. C           IOP    - ONE OF THE VALUES -3,-2,-1,0,1,2,3 GIVEN AS INPUT
  50. C                    AND USED AS DECISION PARAMETER.
  51. C           EPS    - SINGLE PRECISION INPUT VALUE USED AS RELATIVE
  52. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANT DIGITS.
  53. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  54. C                     IER=0  - NO ERROR,
  55. C                     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT
  56. C                              PARAMETERS M,MUD,IOP (SEE REMARKS),
  57. C                              OR BECAUSE OF A NONPOSITIVE RADICAND AT
  58. C                              SOME FACTORIZATION STEP,
  59. C                              OR BECAUSE OF A ZERO DIAGONAL ELEMENT
  60. C                              AT SOME DIVISION STEP.
  61. C                     IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  62. C                              CANCE INDICATED AT FACTORIZATION STEP K+1
  63. C                              WHERE RADICAND WAS NO LONGER GREATER
  64. C                              THAN EPS*A(K+1,K+1).
  65. C
  66. C        REMARKS
  67. C           UPPER PART OF SYMMETRIC BAND MATRIX A CONSISTING OF MAIN
  68. C           DIAGONAL AND MUD UPPER CODIAGONALS (RESP. BAND MATRIX TU
  69. C           CONSISTING OF MAIN DIAGONAL AND MUD UPPER CODIAGONALS)
  70. C           IS ASSUMED TO BE STORED IN COMPRESSED FORM, I.E. ROWWISE
  71. C           IN TOTALLY NEEDED M+MUD*(2M-MUD-1)/2 SUCCESSIVE STORAGE
  72. C           LOCATIONS. ON RETURN UPPER BAND FACTOR TU (ON THE LOCATIONS
  73. C           OF A) IS STORED IN THE SAME WAY.
  74. C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
  75. C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN RESULT MATRIX
  76. C           INVERSE(A)*R OR INVERSE(TU)*R OR INVERSE(TRANSPOSE(TU))*R
  77. C           IS STORED COLUMNWISE TOO ON THE LOCATIONS OF R.
  78. C           INPUT PARAMETERS M, MUD, IOP SHOULD SATISFY THE FOLLOWING
  79. C           RESTRICTIONS     MUD NOT LESS THAN ZERO,
  80. C                            1+MUD NOT GREATER THAN M,
  81. C                            ABS(IOP) NOT GREATER THAN 3.
  82. C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
  83. C           RESTRICTIONS ARE NOT SATISFIED.
  84. C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
  85. C           PARAMETERS ARE SATISFIED, IF RADICANDS AT ALL FACTORIZATION
  86. C           STEPS ARE POSITIVE AND/OR IF ALL DIAGONAL ELEMENTS OF
  87. C           UPPER BAND FACTOR TU ARE NONZERO.
  88. C
  89. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  90. C           NONE
  91. C
  92. C        METHOD
  93. C           FACTORIZATION IS DONE USING CHOLESKY-S SQUARE-ROOT METHOD,
  94. C           WHICH GENERATES THE UPPER BAND MATRIX TU SUCH THAT
  95. C           TRANSPOSE(TU)*TU=A. TU IS RETURNED AS RESULT ON THE
  96. C           LOCATIONS OF A. FURTHER, DEPENDENT ON THE ACTUAL VALUE OF
  97. C           IOP, DIVISION OF R BY TRANSPOSE(TU) AND/OR TU IS PERFORMED
  98. C           AND THE RESULT IS RETURNED ON THE LOCATIONS OF R.
  99. C           FOR REFERENCE, SEE H. RUTISHAUSER, ALGORITHMUS 1 - LINEARES
  100. C           GLEICHUNGSSYSTEM MIT SYMMETRISCHER POSITIV-DEFINITER
  101. C           BANDMATRIX NACH CHOLESKY - , COMPUTING (ARCHIVES FOR
  102. C           ELECTRONIC COMPUTING), VOL.1, ISS.1 (1966), PP.77-78.
  103. C
  104. C     ..................................................................
  105. C
  106.       SUBROUTINE DMCHB(R,A,M,N,MUD,IOP,EPS,IER)
  107. C
  108. C
  109.       DIMENSION R(1),A(1)
  110.       DOUBLE PRECISION TOL,SUM,PIV,R,A
  111. C
  112. C        TEST ON WRONG INPUT PARAMETERS
  113.       IF(IABS(IOP)-3)1,1,43
  114.     1 IF(MUD)43,2,2
  115.     2 MC=MUD+1
  116.       IF(M-MC)43,3,3
  117.     3 MR=M-MUD
  118.       IER=0
  119. C
  120. C        MC IS THE MAXIMUM NUMBER OF ELEMENTS IN THE ROWS OF ARRAY A
  121. C        MR IS THE INDEX OF THE LAST ROW IN ARRAY A WITH MC ELEMENTS
  122. C
  123. C     ******************************************************************
  124. C
  125. C        START FACTORIZATION OF MATRIX A
  126.       IF(IOP)24,4,4
  127.     4 IEND=0
  128.       LLDST=MUD
  129.       DO 23 K=1,M
  130.       IST=IEND+1
  131.       IEND=IST+MUD
  132.       J=K-MR
  133.       IF(J)6,6,5
  134.     5 IEND=IEND-J
  135.     6 IF(J-1)8,8,7
  136.     7 LLDST=LLDST-1
  137.     8 LMAX=MUD
  138.       J=MC-K
  139.       IF(J)10,10,9
  140.     9 LMAX=LMAX-J
  141.    10 ID=0
  142.       TOL=A(IST)*EPS
  143. C
  144. C        START FACTORIZATION-LOOP OVER K-TH ROW
  145.       DO 23 I=IST,IEND
  146.       SUM=0.D0
  147.       IF(LMAX)14,14,11
  148. C
  149. C        PREPARE INNER LOOP
  150.    11 LL=IST
  151.       LLD=LLDST
  152. C
  153. C        START INNER LOOP
  154.       DO 13 L=1,LMAX
  155.       LL=LL-LLD
  156.       LLL=LL+ID
  157.       SUM=SUM+A(LL)*A(LLL)
  158.       IF(LLD-MUD)12,13,13
  159.    12 LLD=LLD+1
  160.    13 CONTINUE
  161. C        END OF INNER LOOP
  162. C
  163. C        TRANSFORM ELEMENT A(I)
  164.    14 SUM=A(I)-SUM
  165.       IF(I-IST)15,15,20
  166. C
  167. C        A(I) IS DIAGONAL ELEMENT. ERROR TEST.
  168.    15 IF(SUM)43,43,16
  169. C
  170. C        TEST ON LOSS OF SIGNIFICANT DIGITS AND WARNING
  171.    16 IF(SUM-TOL)17,17,19
  172.    17 IF(IER)18,18,19
  173.    18 IER=K-1
  174. C
  175. C        COMPUTATION OF PIVOT ELEMENT
  176.    19 PIV=DSQRT(SUM)
  177.       A(I)=PIV
  178.       PIV=1.D0/PIV
  179.       GO TO 21
  180. C
  181. C        A(I) IS NOT DIAGONAL ELEMENT
  182.    20 A(I)=SUM*PIV
  183. C
  184. C        UPDATE ID AND LMAX
  185.    21 ID=ID+1
  186.       IF(ID-J)23,23,22
  187.    22 LMAX=LMAX-1
  188.    23 CONTINUE
  189. C
  190. C        END OF FACTORIZATION-LOOP OVER K-TH ROW
  191. C        END OF FACTORIZATION OF MATRIX A
  192. C
  193. C     ******************************************************************
  194. C
  195. C        PREPARE MATRIX DIVISIONS
  196.       IF(IOP)24,44,24
  197.    24 ID=N*M
  198.       IEND=IABS(IOP)-2
  199.       IF(IEND)25,35,25
  200. C
  201. C     ******************************************************************
  202. C
  203. C        START DIVISION BY TRANSPOSE OF MATRIX TU (TU IS STORED IN
  204. C        LOCATIONS OF A)
  205.    25 IST=1
  206.       LMAX=0
  207.       J=-MR
  208.       LLDST=MUD
  209.       DO 34 K=1,M
  210.       PIV=A(IST)
  211.       IF(PIV)26,43,26
  212.    26 PIV=1.D0/PIV
  213. C
  214. C        STA-T BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  215.       DO 30 I=K,ID,M
  216.       SUM=0.D0
  217.       IF(LMAX)30,30,27
  218. C
  219. C        PREPARE INNER LOOP
  220.    27 LL=IST
  221.       LLL=I
  222.       LLD=LLDST
  223. C
  224. C        START INNER LOOP
  225.       DO 29 L=1,LMAX
  226.       LL=LL-LLD
  227.       LLL=LLL-1
  228.       SUM=SUM+A(LL)*R(LLL)
  229.       IF(LLD-MUD)28,29,29
  230.    28 LLD=LLD+1
  231.    29 CONTINUE
  232. C        END OF INNER LOOP
  233. C
  234. C        TRANSFORM ELEMENT R(I)
  235.    30 R(I)=PIV*(R(I)-SUM)
  236. C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  237. C
  238. C        UPDATE PARAMETERS LMAX, IST AND LLDST
  239.       IF(MC-K)32,32,31
  240.    31 LMAX=K
  241.    32 IST=IST+MC
  242.       J=J+1
  243.       IF(J)34,34,33
  244.    33 IST=IST-J
  245.       LLDST=LLDST-1
  246.    34 CONTINUE
  247. C
  248. C        END OF DIVISION BY TRANSPOSE OF MATRIX TU
  249. C
  250. C     ******************************************************************
  251. C
  252. C        START DIVISION BY MATRIX TU (TU IS STORED ON LOCATIONS OF A)
  253.       IF(IEND)35,35,44
  254.    35 IST=M+(MUD*(M+M-MC))/2+1
  255.       LMAX=0
  256.       K=M
  257.    36 IEND=IST-1
  258.       IST=IEND-LMAX
  259.       PIV=A(IST)
  260.       IF(PIV)37,43,37
  261.    37 PIV=1.D0/PIV
  262.       L=IST+1
  263. C
  264. C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  265.       DO 40 I=K,ID,M
  266.       SUM=0.D0
  267.       IF(LMAX)40,40,38
  268.    38 LLL=I
  269. C
  270. C        START INNER LOOP
  271.       DO 39 LL=L,IEND
  272.       LLL=LLL+1
  273.    39 SUM=SUM+A(LL)*R(LLL)
  274. C        END OF INNER LOOP
  275. C
  276. C        TRANSFORM ELEMENT R(I)
  277.    40 R(I)=PIV*(R(I)-SUM)
  278. C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  279. C
  280. C        UPDATE PARAMETERS LMAX AND K
  281.       IF(K-MR)42,42,41
  282.    41 LMAX=LMAX+1
  283.    42 K=K-1
  284.       IF(K)44,44,36
  285. C
  286. C        END OF DIVISION BY MATRIX TU
  287. C
  288. C     ******************************************************************
  289. C
  290. C        ERROR EXIT IN CASE OF WRONG INPUT PARAMETERS OR PIVOT ELEMENT
  291. C        LESS THAN OR EQUAL TO ZERO
  292.    43 IER=-1
  293.    44 RETURN
  294.       END
  295.